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When common factors strongly influence two power-law cross-correlated time series recorded in 
complex natural or social systems, using classic detrended cross-correlation analysis (DCCA) without 
considering these common factors will bias the results. We use detrended partial cross-correlation 
analysis (DPXA) to uncover the intrinsic power-law cross-correlations between two simultaneously 
recorded time series in the presence of nonstationarity after removing the effects of other time series 
acting as common forces. The DPXA method is a generalization of the detrended cross-correlation 
analysis that takes into account partial correlation analysis. We demonstrate the method by using 
bivariate fractional Brownian motions contaminated with a fractional Brownian motion. We find 
that the DPXA is able to recover the analytical cross Hurst indices, and thus the multi-scale DPXA 
coefficients are a viable alternative to the conventional cross-correlation coefficient. We demonstrate 
the advantage of the DPXA coefficients over the DCCA coefficients by analyzing contaminated 
bivariate fractional Brownian motions. We calculate the DPXA coefficients and use them to extract 
the intrinsic cross-correlation between crude oil and gold futures by taking into consideration the 
impact of the US dollar index. We develop the multifractal DPXA (MF-DPXA) method in order 
to generalize the DPXA method and investigate multifractal time series. We analyze multifractal 
binomial measures masked with strong white noises and find that the MF-DPXA method quantifies 
the hidden multifractal nature while the MF-DCCA method fails. 

PACS numbers: 89.75.Da, 05.45.Tp, 05.45.Df, 05.40.-a 


I. INTRODUCTION 

Complex systems with interacting constituents are 
ubiquitous in nature and society. To understand the 
microscopic mechanisms of emerging statistical laws of 
complex systems, one records and analyzes time series 
of observable quantities. These time series are usually 
nonstationary and possess long-range power-law cross¬ 
correlations. Examples include the velocity, temperature, 
and concentration fields of turbulent flows embedded in 
the same space as joint multifractal measures [ 1 U 4 I , to¬ 
pographic indices and crop yield in agronomy [3|, 1^ , tem¬ 
poral and spatial seismic data nitrogen dioxide and 
ground-level ozone @ , heart rate variability and brain ac¬ 
tivity in healthy humans Q, sunspot numbers and river 
flow fluctuations Q, wind patterns and land surface air 
temperatures @ , traffic flows m and traffic signals [HI: 
self-affine time series of taxi accidents [l^, and econo¬ 
physical variables (l3l - [r^ . 

A variety of methods have been used to investigate 
the long-range power-law cross-correlations between two 
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nonstationary time series. The earliest was joint multi- 
fractal analysis to study the cross-multifractal nature of 
two joint multifractal measures through the scaling be¬ 
haviors of the joint moments [H, H, [l 8 l - [^ . which is a 
multifractal cross-correlation analysis based on the par¬ 
tition function approach (MF-X-PF) [2l|. Over the past 
decade, detrended cross-correlation analysis (DCCA) has 
become the most popular method of investigating the 
long-range power-law cross correlations between two non¬ 
stationary time series H, [ 22 l - [^ , and this method has 
numerous variants [25l - l33l |. Statistical tests can be used 
to measure these cross correlations There is 

also a group of multifractal detrended fluctuation analy¬ 
sis (MF-DCCA) methods of analyzing multifractal time 
series, e.g ., MF-X-DFA [13, MF-X-DMA [H, and MF- 
HXA da . 

The observed long-range power-law cross-correlations 
between two time series may not be caused by their in¬ 
trinsic relationship but by a common third driving force 
or by common external factors [dol - ld^ . If the influence 
of the common external factors on the two time series 
are additive, we can use partial correlation to measure 
their intrinsic relationship [43j |. To extract the intrin¬ 
sic long-range power-law cross-correlations between two 
time series affected by common driving driving forces, we 
previously developed and used detrended partial cross- 
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correlation analysis (DPXA) and studied the DPXA ex¬ 
ponents of variable cases, combining the ideas of de¬ 
trended cross-correlation analysis and partial correlation 
(ill . In Ref. the DPXA method has been proposed 
independently, focussing on the DPXA coefficient. 

Here we provide a general framework for the DPXA 
and MF-DPXA methods that is applicable to various ex¬ 
tensions, including different detrending approaches and 
higher dimensions. We adopt two well-established math¬ 
ematical models (bivariate fractional Brownian motions 
and multifractal binomial measures) in our numerical ex¬ 
periments, which have known analytical expressions, and 
demonstrate how the (MF-)DPXA methods is superior to 
the corresponding (MF-)DCCA methods. 


II. DETRENDED PARTIAL 
CROSS-CORRELATION ANALYSIS 


A. DPXA exponent 


Consider two stationary time series {x(t) : t = 

I,-- - ,T} and {y{t) : t = I,-- - ,T} that depend on a 
sequence of time series {zi{t) : t = 1,2, ••• ,T} with 
i = 1, • ■ • , n. Each time series is covered with Mg = \T/s] 
non-overlapping windows of size s. Consider the rith box 
[/« -I- 1, + s], where ly = (u — l)s. We calibrate the two 

linear regression models for Xy and yy respectively, 

J Xy - ZyfBx^V “t“ Xx^y /. \ 

( Vv — '^vf^y,v “t” ^y,v 


where = [xi^+i,.. .,xi^+s]'^, y„ = [yi,+i, ■ ■. , 2 /;„+s]'^, 
Vx^y and Vy^y are the vectors of the error term, and 


Zy — 



^ Zl{ly + 1 ) 

\ Zl{ly + S) 


Zpi^ly -f 1) 
Zpi^ly A 5) 


( 2 ) 


is the matrix of the p external forces in the uth box, 
where is the transform of x. Equation © gives the 
estimates /3x,v and /3y^y of the p-dimensional parameter 
vectors f3x,y and f3y^y and the sequence of error terms. 


j ^X,V - 

L ^y,v — yv '^vf3y,v 


We obtain the disturbance profiles, i.e., 


Rx,v{k) = rx{ly+j) 

Ry,v{k) = = l '^y{^v + j) 


(3) 

(4) 


where k = 1 , • • • , s. 

We assume that the local trend functions of Rx^y and 
Rx^y are Rx,v and Ry^y, respectively. The detrended par¬ 
tial cross-correlation in each window is then calculated. 




Rx,v{k') 


R’Vj'xi.k') Ry^xik') , 


( 5 ) 


and the second-order detrended partial cross-correlation 
is calculated, 


Fxy.z{2,s) 


1 


TO — 1 


m 


1/2 




( 6 ) 


If there are intrinsic long-range power-law cross¬ 
correlations between x and y, we expect the scaling rela¬ 
tion. 


Fxy.z{‘2.,s) 


(7) 


There are many ways of determining Rx^v and Ry^y. 
The local detrending functions could be polynomials 
4^ Id^ . moving averages (4§ - l^ . or other possibilities 
52l| . To distinguish the different detrending methods, we 
label the corresponding DPXA variants as, e.g., PX-DFA 
and PX-DMA. When the moving average is used as the 
local detrending function, the window size of the moving 
averages must be the same as the covering window size s 


To measure the validity of the DPXA method, we per¬ 
form numerical experiments using an additive model for 
x and y, i.e.. 


J x{t) = l3x,o + Pxz{t) + rx{t) 

\ y{t) = Pyfi + Pyzit) + ry{t) ’ 

where z{t) is a fractional Gaussian noise with Hurst index 
Hz, and Xx and Xy are the incremental series of the two 
components of a bivariate fractional Brownian motion 
(BEBMs) with Hurst indices Hy^ and Hy^ The 

properties of multivariate fractional Brownian motions 
have been extensively studied In particular, it 

has been proven that the Hurst index Hy^y^ of the cross¬ 
correlation between the two components is 

Hy^y^ = {Hy^+Hy^)/2. (9) 

This property allows us to assess how the proposed 
method perform. We can obtain the hxy of x and y using 
the DCCA method and the hxy.z of Xx and Xy using the 
DPXA method. Our numerical experiments show that 
Hy^cTy = ky^yy = kxy.z ^ hxy Wo USB H for theoretical 
or true values and h for numerical estimates. 

In the simulations we set j3xfl = 2, Px = 3, j3yfi = 2, 
and Py = 3 in the model based on Eq. ([5]). Three 
Hurst indices Hy ^, Hy ^, and are input arguments 
and vary from 0.1 to 0.95 at 0.05 intervals. Because 
Xx and Xy are symmetric, we set Hy^ < Hy^, result¬ 
ing in X 18 = 3078 triplets of {Hy^, Hy^, Hz). 

The BFBMs are simulated using the method described 
in Ref. [5^ [H^, and the FBMs are generated using a 
rapid wavelet-based approach [s^]- The length of each 
time series is 65536. For each {Hy^, Hy^, Hz) triplet we 
conduct 100 simulations. We obtain the Hurst indices 
for the simulated time series ^ ^ z, x, and y using de¬ 
trended fluctuation analysis • The average values 
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FIG. 1. (color online.) Detrended partial cross-correlation exponents, (a-b) Dependence of the average Hurst indices hr^ and 
hr^ of the two components of the generated bivariate fractional Brownian motions on the input Hurst indices Hr^ and Hr ^. (c) 
Dependence of the average Hurst index of the generated univariate fractional Brownian motions on the input Hurst index 
Hz. (d) Dependence of hx on Hz for different hr^ values, (e) Relative error Ahr^ry = hr^ry — Hr^ry (f) Relative error Ahxy.z 
between the DPXA estimate {hxy.z)z and the true value hr^ry as a function of hr^ and hry 


hr , hz, hx, and hy over 100 realizations are cal¬ 
culated for further analysis, which are shown in Fig. [T] 
A linear regression between the output and input Hurst 
indices in Fig. [Ha-c) yields {hr^) = 0.009 + O.OOOHr^j 
{hry) = 0.009 + Q.mQHry, and {K) = 0.010 + 0.991iJ^, 
suggesting that the generated FBMs have Hurst indices 
equal to the input Hurst indices. Figure [ijd) shows 
that when hr,, < hz, hx is close to hz. When it is not, 
hz < hx < hr„ ■ 

Figure [He) shows that hr y,ry = {hr,, +hry)/2. Because 
hr„ ss Hr^ and Ky ^ Hry [see Fig. |Il[a)-(b)], we verify 
numerically that 


hr„ry^Hr„ry- ( 10 ) 

Note also that hxy ~ {hx + hy)/2, and that hxy.z is a 
function of hr „, hry and hz ■ A simple linear regression 
gives 


hxy.z = 0.003 + 0.509h^, -b 0.493h^^ + 0.012h^, (11) 

which indicates that the DPXA method can be used to 
extract the intrinsic cross-correlations between the two 
time series x and y when they are influenced by a common 
factor z. We calculate the average {hxy.z)z over different 
Hz: and then find the relative error 


Ahxy.z — 


{hxy.z) z hr, 


( 12 ) 


Figure [TJf) shows the results for different combinations 
of hry, and hry. Although in most cases we see that 
Ahxy.z ^ 0.05, when both hr^ and hr approach 0, 
Ahxy.z increases. When hr^ = hry = 0.11, Ahxy.z = 
0.192, and when hr^ = 0.11 and hry = 0.16, Ahxy.z = 
0.113. For all other points of {hry,hry), the relative er¬ 
rors Ahxy.z are less than 0.10. 


B. DPXA coefficient 


In a way similar to detrended cross-correlation coef¬ 
ficients 1^, we define the detrended partial cross¬ 
correlation coefficient (or DPXA coefficient) as 


PDPXa(s) = Pxy.z{s) 


^xy.z{‘^^ '®) 

Fx:z{2,s)Fy,z{2,s)' 


(13) 


As in the DCCA coefficient (sT], we also find — 1 < 
Pdpxa(s) < 1 for DPXA. The DPXA coefficient indi¬ 
cates the intrinsic cross-correlations between two non¬ 
stationary series. 

We use the mathematical model in Eq. (|8]) with the co¬ 
efficients /3x,o = /?y,o = 2 and /3x,i = Py,i = 3 to demon¬ 
strate how the DPXA coefficient outperforms the DCCA 
coefficient. The two components Tx and Vy of the BFBM 
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FIG. 2. (color online.) Detrended partial cross-correlation 
coefficients, (a) Performance of different methods by compar¬ 
ing three cross-correlation coefficients Px,y, Prx,ry and px,y.z 
of the mathematical model in Eq. ([8|. (b) Estimation and 
comparison of the cross-correlation levels between the two re¬ 
turn time series (o) and two volatility time series (□) of crude 
oil and gold when including and excluding the influence of the 
USD index. 


have very small Hurst indices Hr^ = 0.1 = Hr^ = 0.1 
and their correlation coefficient is p = 0.7, and the driv¬ 
ing FBM force z has a large Hurst index = 0.95. 
Figure UKa) shows the resulting cross-correlation coeffi¬ 
cients at different scales. The DCCA coefficients Pra,,ry 
between the generated and Vy time series overestimate 
the true value p = 0.7. Because the influence of z on 
and ry is very strong, the behaviors of x and y are dom¬ 
inated by z, and the cross-correlation coefficient Px,y{s) 
is close to 1 when s is small and approaches 1 when s 
us large. In contrast, the DPXA coefficients Px,y:z are 
in good agreement with the true value p = 0.7. Note 
that the DPXA method better estimates Tx and ry than 
the DCCA method, since the Pry,,ry curve deviates more 
from the horizontal line p = 0.7 than the Px,y:z curve, 
especially at large scales. 


To illustrate the method with an example from finance, 
we use it to estimate the intrinsic cross-correlation levels 
between the futures returns and the volatilities of crude 
oil and gold. It is well-documented that the returns of 
crude oil and gold futures are correlated [s^, and that 
both commodities are influenced by the USD index [^ . 
The data samples contain the daily closing prices of gold, 
crude oil, and the USD index from 4 October 1985 to 31 
October 2012. Figure [2])b) shows that both the DCCA 
and DPXA coefficients of returns exhibit an increasing 
trend with respect to the scale s, and that the two types 
of coefficient for the volatilities do not exhibit any evident 
trend. For both financial variables. Fig. [2Db) shows that 

Pg.o:d(s) < Pg,o(s) (14) 

for different scales. Although this is similar to the 
result between ordinary partial correlations and cross¬ 
correlations , the DPXA coefficients contain more in¬ 
formation than the ordinary partial correlations since the 
former indicate the partial correlations at multiple scales. 


III. MULTIFRACTAL DETRENDED PARTIAL 
CROSS-CORRELATION ANALYSIS 


An extension of the DPXA for multifractal time series, 
notated MF-DPXA, can be easily implemented. When 
MF-DPXA is implemented with DFA or DMA, we no¬ 
tate it MF-PX-DFA or MF-PX-DMA. The gth order de¬ 
trended partial cross-correlation is calculated 


TO — 1 


Fxy.z (Qj s) — 
when g 0, and 

Fxy.z{^,s) = exp 


1 m 


-[ 1/9 


■^ln|F„(s)| 


V — 1 


We then expect the scaling relation 

Fxy.z{q,s) ^ . 


(15) 


(16) 


(17) 


According to the standard multifractal formalism, the 
multifractal mass exponent T{q) can be used to charac¬ 
terize the multifractal nature, i.e.. 


’^xy.zisi) — Q^xy.z^q) F f (1^) 

where Df is the fractal dimension of the geometric sup¬ 
port of the multifractal measure . We use Df = 1 for 
our time series analysis. If the mass exponent r(g) is a 
nonlinear function of q, the signal is multifractal. We use 
the Legendre transform to obtain the singularity str eng th 
function a{q) and the multifractal spectrum f{a) [63[ 

f OCf^yM) = dT:,y..f,{q)/dq 
\ fxy:z{q) — qC^xy.z '^xy:z{q) 

























5 



FIG. 3. (color online.) Multifractal detrended partial cross-correlation analysis of two binomial measures contaminated by 
Gaussian noise with very low signal-to-noise ratio, (a-c) The segments of the binomial signal rx{t) with = 0.3, the Gaussian 
noise z{t), and the “observed” signal x{t). (d-f) Power-law dependence of the fluctuations Fxy{q,s), Fxy:z{q, s), Fr^ryiq, s) on 
the scale s for q = —4, 0 and 4 (top down). The values of F^y.z and Frxry have been multiplied by 10®. (g) Multifractal mass 
exponents rxy{q), Txy.z{q) and rrxry{q), with the theoretical curve Tr^ry shown as a continuous line, (d) Multifractal spectra 
fxy(a}, fxy.z{a), fr^ryia) of the singularity strength oi. The continuous curve is the theoretical spectrum Fr^ryia)- 


To test the performance of MF-DPXA, we construct 
two binomial measures {rx{t) : t = 1 , 2 ,-•• , 2 ^} and 
{ry{t) : t = 1, 2 , • • ■ , 2 ^'} from then-model with known 
analytic multifractal properties j64| . and contaminate 
them with Gaussian noise. We generate the binomial 
measure iteratively [s^ by using the multiplicative fac¬ 
tors Px = 0.3 for Tx and Py = 0.4 for Vy. The contami¬ 
nated signals are a: = 2 -|- 3z -I- and y = 2 + Zz + Vy. 
Figures [Sla)-(c) show that the signal-to-noise ratio is 
of order 0(10“®). Figures |3Id)-(f) show a power-law 
dependence between the fluctuation functions and the 
scale, in which it is hard to distinguish the three curves 
of Fxy Figure E^g) shows that for x{t) and y{t), the 
Txy (<?) function an approximate straight line and that the 
corresponding fxyio:) spectrum is very narrow and con¬ 
centrated around a = 0.5. These observations are trivial 
because x{t) and y(t) are Gaussian noise with the Hurst 
indices Hx = Hy = 0.5, and the multifractal detrended 
cross-correlation analysis fails to uncover any mul- 
tifractality. On the contrary, we find that Txy.ziq) ~ 
Fxryiq) - TrxVyiq) and fxy:z{a) - fr^Vyia) ~ 

Thus the MF-DPXA method successfully reveals the in¬ 
trinsic multifractal nature between rx{t) and ry(t) hidden 
in x{t) and y{t). 


IV. SUMMARY 

In summary, we have studied the performances of 
DPXA exponents, DPXA coefficients, and MF-DPXA 
using bivariate fractional Brownian motions contami¬ 
nated by a fractional Brownian motion and multifractal 
binomial measures contaminated by white noise. These 
mathematical models are appropriate here because their 
analytical expressions are known. We have demonstrated 
that the DPXA methods are capable of extracting the 
intrinsic cross-correlations between two time series when 
they are influenced by common factors, while the DGGA 
methods fail. 

The methods discussed are intended for multivariate 
time series analysis, but they can also be generalized to 
higher dimensions [^,[53,[^|^. We can also use lagged 
cross-correlations in these methods [13, H^. Although 
comparing the performances of different methods is al¬ 
ways important [b^, different variants of a method can 
produce different outcomes when applied to different sys¬ 
tems. For instance, one variant that outperforms other 
variants under the setting of certain stochastic processes 
is not necessary the best performing method for other 
systems [s^. We argue that there are still a lot of open 
questions for the big family of DFA, DMA, DCCA and 
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